//
//  KinematicsCalculator.m
//  MacHeliosSim
//
//  Created by Paul Thompson on 8/12/10.
//  Copyright 2010 __MyCompanyName__. All rights reserved.
//

#import "KinematicsCalculator.h"

@implementation KinematicsCalculator
-(id)initWithReaction:(Reaction *)theReaction
{
	myNuclearReaction = theReaction;
	stuffToReturn = [[KinCalcReturnData alloc] init];
	return self;	
}

-(KinCalcReturnData *)labEjectileEnergyForAngle:(float)theAngle
{
	float centerMassTotalEnergy;
	float centerMassEjectileEnergy;
	float kEjectile;	//	Constant defined as centerMassFrameVelocity/centerMassParticle3Velocity
						//	particle 3 is the ejectile, particle 4 is the product, e.g., proton and 
						//	20O respectively for the test case.
	
	float labEjectileEnergy;
	float labBeamEnergy;
	float massBeam, massTarget, massEjectile, massProduct;
	float massConverter = 931.494028;  // REMEMBER TO DIVIDE ALL DELTAM BY 1000!!!!!!!
	float qValue;
	
	
	// keep masses in MeV.  Note that all mass defects (deltaM) are in KeV, not MeV, and thus must be divided by 1000....
	// that's really stupid.  Why don't I just change them all?
	massBeam = myNuclearReaction.beamNuclide.nucleonCount*massConverter+myNuclearReaction.beamNuclide.deltaM/1000.0;
	massTarget = myNuclearReaction.targetNuclide.nucleonCount*massConverter+myNuclearReaction.targetNuclide.deltaM/1000.0;
	massEjectile = myNuclearReaction.ejectileNuclide.nucleonCount*massConverter+myNuclearReaction.ejectileNuclide.deltaM/1000.0;
	massProduct = myNuclearReaction.productNuclide.nucleonCount*massConverter+myNuclearReaction.productNuclide.deltaM/1000.0;
	qValue = myNuclearReaction.qValue;
	labBeamEnergy = myNuclearReaction.beamEnergy; //Also in MeV
	
	centerMassTotalEnergy = labBeamEnergy*massTarget/(massBeam+massTarget);

	kEjectile = sqrt( (massBeam*massEjectile*centerMassTotalEnergy*(massEjectile+massProduct))
					 /(massTarget*massProduct*(centerMassTotalEnergy+qValue)*(massBeam+massTarget)) );//ejectile = particle 3


	centerMassEjectileEnergy = (massProduct*(centerMassTotalEnergy+qValue))
							  /(massEjectile + massProduct);
											
	labEjectileEnergy = centerMassEjectileEnergy*(1.0 + pow(kEjectile, 2.0) + 2.0 * kEjectile * cos(theAngle));
		
		
	[stuffToReturn setEnergy:labEjectileEnergy];
	
	float tempLabAngle = (atan(sin(theAngle)/(cos(theAngle)+kEjectile)));
	if( signCompare(tempLabAngle, theAngle) )
	{
		[stuffToReturn setLabAngle:tempLabAngle];
	} else if (theAngle > 0.0) {
		[stuffToReturn setLabAngle:(pi + tempLabAngle)];
	} else { // theAngle is negative
		[stuffToReturn setLabAngle:(-pi + tempLabAngle)]; 
	}

	
	if (kEjectile < 1.0) {
		[stuffToReturn setMaxLabAngle:pi];
	} else if (kEjectile == 1.0) {
		[stuffToReturn setMaxLabAngle:pi/2.0];
	} else {  // kEjectile > 1.0
		[stuffToReturn setMaxLabAngle:(atan(1.0/(sqrt(pow(kEjectile, 2.0)-1.0))))];
	}
	
	[stuffToReturn setCmAngle:theAngle];
	
#ifdef DEBUG
	NSLog(@"Ejectile Energy for cm angle %f degrees, lab angle %f degrees, is %f MeV",[stuffToReturn cmAngleInDegrees], [stuffToReturn labAngleInDegrees], labEjectileEnergy);				 
#endif
	return stuffToReturn;
}


-(KinCalcReturnData *)labProductEnergyForAngle:(float)theAngle
{
	float centerMassTotalEnergy;
	float centerMassProductEnergy;
	float kProduct;	//	Constant defined as centerMassFrameVelocity/centerMassParticle4Velocity
					//	particle 3 is the ejectile, particle 4 is the product, e.g., proton and
					//	20O respectively for the test case.
	
	float labProductEnergy;
	float labBeamEnergy;
	float massBeam, massTarget, massEjectile, massProduct;
	float massConverter = 931.494028;  // REMEMBER TO DIVIDE ALL DELTAM BY 1000!!!!!!!
	float qValue;
	
	// keep masses in MeV.  Note that all mass defects (deltaM) are in KeV, not MeV, and thus must be divided by 1000....
	// that's really stupid.  Why don't I just change them all?
	
	massBeam = myNuclearReaction.beamNuclide.nucleonCount*massConverter+myNuclearReaction.beamNuclide.deltaM/1000.0;
	massTarget = myNuclearReaction.targetNuclide.nucleonCount*massConverter+myNuclearReaction.targetNuclide.deltaM/1000.0;
	massEjectile = myNuclearReaction.ejectileNuclide.nucleonCount*massConverter+myNuclearReaction.ejectileNuclide.deltaM/1000.0;
	massProduct = myNuclearReaction.productNuclide.nucleonCount*massConverter+myNuclearReaction.productNuclide.deltaM/1000.0;
	qValue = myNuclearReaction.qValue;
	labBeamEnergy = myNuclearReaction.beamEnergy; //Also in MeV
	
	centerMassTotalEnergy = labBeamEnergy*massTarget/(massBeam+massTarget);


	kProduct = sqrt( (massBeam*massProduct*centerMassTotalEnergy*(massEjectile+massProduct))
					/(massTarget*massEjectile*(centerMassTotalEnergy+qValue)*(massBeam+massTarget)) );//product = particle 4
	
	centerMassProductEnergy =  (massEjectile*(centerMassTotalEnergy+qValue))
	/(massEjectile+massProduct);
	
	labProductEnergy = centerMassProductEnergy * (1.0 + pow(kProduct, 2.0) + 2.0 * kProduct * cos(theAngle));

	
	
	[stuffToReturn setEnergy:labProductEnergy];

	float tempLabAngle = (atan(sin(theAngle)/(cos(theAngle)+kProduct)));
	if( signCompare(tempLabAngle, theAngle) )
	{
		[stuffToReturn setLabAngle:tempLabAngle];
	} else if (theAngle > 0.0) {
		[stuffToReturn setLabAngle:(pi + tempLabAngle)];
	} else { // theAngle is negative
		[stuffToReturn setLabAngle:(-pi + tempLabAngle)]; 
	}
	
	if (kProduct < 1.0) {
		[stuffToReturn setMaxLabAngle:pi];
	} else if (kProduct == 1.0) {
		[stuffToReturn setMaxLabAngle:pi/2.0];
	} else {  // kProduct > 1.0
		[stuffToReturn setMaxLabAngle:(atan(1.0/(sqrt(pow(kProduct, 2.0)-1.0))))];
	}
	
	[stuffToReturn setCmAngle:theAngle];
	 
	
	
#ifdef DEBUG
	NSLog(@"Product Energy for cm angle %f degrees, lab angle %f degrees, is %f MeV",[stuffToReturn cmAngleInDegrees], [stuffToReturn labAngleInDegrees], labProductEnergy);
#endif
	return stuffToReturn;
	
}
								
	
											   

@end


bool signCompare (float someFloat, float anotherFloat)
{
	return	( 
			 (someFloat < 0 && anotherFloat < 0) 
			 || (someFloat > 0 && anotherFloat > 0)
			 || (someFloat == 0 && anotherFloat == 0)
			 );
}
